﻿using System;
using System.Collections.Generic;

namespace Pro2
{
    enum TransOption
    {
        Three = 1,  // 三度带
        Six = 2     // 六度带
    }

    class CoordTrans
    {
        public Ellipsoid Ellipsoid { set; get; }

        /// <summary>
        /// 默认构造函数，使用克拉索夫斯基椭球
        /// </summary>
        public CoordTrans()
        {
            Ellipsoid = new Ellipsoid(EllipsoidType.KRASOVSKI);
        }

        public CoordTrans(EllipsoidType type)
        {
            Ellipsoid = new Ellipsoid(type);
        }

        /// <summary>
        /// 高斯投影正算
        /// </summary>
        /// <param name="B"></param>
        /// <param name="L"></param>
        /// <param name="L0"></param>
        /// <param name="n"></param>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="Y"></param>
        public void BLtoXY(double B, double L, double L0, int n, out double x, out double y, out double Y)
        {
            // 计算子午线弧长
            double[] alphaCoeff = new double[6];
            ComputeAlphaCoeff(ref alphaCoeff);
            double X_zwx = alphaCoeff[0] * B + alphaCoeff[1] * Math.Sin(2 * B) + alphaCoeff[2] * Math.Sin(4 * B) + alphaCoeff[3] * Math.Sin(6 * B) + alphaCoeff[4] * Math.Sin(8 * B) + alphaCoeff[5] * Math.Sin(10 * B);

            // 计算经差
            double l = L - L0;

            // 计算辅助量
            double[] aCoeff = new double[7];
            ComputeACoeff(B, X_zwx, ref aCoeff);

            // 计算自然坐标
            x = aCoeff[0] + aCoeff[2] * Math.Pow(l, 2) + aCoeff[4] * Math.Pow(l, 4) + aCoeff[6] * Math.Pow(l, 6);
            y = aCoeff[1] * l + aCoeff[3] * Math.Pow(l, 3) + aCoeff[5] * Math.Pow(l, 5);

            // 计算通用坐标
            Y = n * 1000000 + y + 500000;
        }

        /// <summary>
        /// 高斯投影反算
        /// </summary>
        /// <param name="X"></param>
        /// <param name="Y"></param>
        /// <param name="option"></param>
        /// <param name="B"></param>
        /// <param name="L"></param>
        public void XYtoBL(double X, double Y, TransOption option,out double B, out double L)
        {
            // 计算中央子午线经度L0
            int n = (int)(Y / 1000000);
            int L0;
            if (option == TransOption.Three)
            {
                L0 = 3 * n;
            }
            else
            {
                L0 = 6 * n - 3;
            }

            // 计算自然坐标（x,y)
            double x = X;
            double y = (int)(Y / 1000000) * 1000000 - 500000;

            // 计算底点纬度
            double[] alphaCoeff = new double[6];
            ComputeAlphaCoeff(ref alphaCoeff);
            double B0 = X / alphaCoeff[0];
            double delta = alphaCoeff[1] * Math.Sin(2 * B0) + alphaCoeff[2] * Math.Sin(4 * B0) + alphaCoeff[3] * Math.Sin(6 * B0) + alphaCoeff[4] * Math.Sin(8 * B0) + alphaCoeff[5] * Math.Sin(10 * B0);
            double Bf = (X - delta) / alphaCoeff[0];
            while (Math.Abs(Bf - B0) < 1e-8)
            {
                B0 = Bf;
                delta = alphaCoeff[1] * Math.Sin(2 * B0) + alphaCoeff[2] * Math.Sin(4 * B0) + alphaCoeff[3] * Math.Sin(6 * B0) + alphaCoeff[4] * Math.Sin(8 * B0) + alphaCoeff[5] * Math.Sin(10 * B0);
                Bf = (X - delta) / alphaCoeff[0];
            }

            // 计算辅助量
            double[] bCoeff = new double[7];
            ComputeBCoeff(Bf, ref bCoeff);

            // 计算（B，L）
            B = bCoeff[0] + bCoeff[2] * Math.Pow(y, 2) + bCoeff[4] * Math.Pow(y, 4) + bCoeff[6] * Math.Pow(y, 6);
            L = bCoeff[1] * y + bCoeff[3] * Math.Pow(y, 3) + bCoeff[5] * Math.Pow(y, 5) + L0;
        }

        /// <summary>
        /// 六度带转三度带
        /// </summary>
        /// <param name="data">输入数据</param>
        /// <returns></returns>
        public Data SixToThree(Data data)
        {
            double B, L, L0, x, y, Y;
            int n;
            XYtoBL(data.X, data.Y_general, TransOption.Six, out B, out L);
            ComputeL0AndN(L, TransOption.Three, out L0, out n);
            BLtoXY(B, L, L0, n, out x, out y, out Y);
            return new Data { ID = data.ID, X = x, Y = y, Y_general = Y };
        }

        /// <summary>
        /// 三度带转六度带
        /// </summary>
        /// <param name="data">输入数据</param>
        /// <returns></returns>
        public Data ThreeToSix(Data data)
        {
            double B, L, L0, x, y, Y;
            int n;
            XYtoBL(data.X, data.Y_general, TransOption.Three, out B, out L);
            ComputeL0AndN(L, TransOption.Six, out L0, out n);
            BLtoXY(B, L, L0, n, out x, out y, out Y);
            return new Data { ID = data.ID, X = x, Y = y, Y_general = Y };
        }

        /// <summary>
        /// 计算邻带坐标
        /// </summary>
        /// <param name="data">输入数据</param>
        /// <param name="option">在三度带还是六度带情况下进行转换</param>
        /// <returns></returns>
        public List<Data> NeighBorCoord(Data data, TransOption option)
        {
            List<Data> datas = new List<Data>();

            double B, L, L0;
            int n;
            XYtoBL(data.X, data.Y_general, option, out B, out L);
            ComputeL0AndN(L, option, out L0, out n);

            // 计算右侧坐标
            double x1, y1, Y1;
            double L1;
            int n1;
            if (option == TransOption.Three)
                L1 = L0 + 3;
            else
                L1 = L0 + 6;
            n1 = n + 1;
            BLtoXY(B, L, L1, n1, out x1, out y1, out Y1);
            datas.Add(new Data { ID = data.ID, X = x1, Y = y1, Y_general = Y1 });

            // 计算左侧坐标
            double x2, y2, Y2;
            double L2;
            int n2;
            if (option == TransOption.Three)
                L2 = L0 - 3;
            else
                L2 = L0 - 6;
            n2 = n - 1;
            BLtoXY(B, L, L2, n2, out x2, out y2, out Y2);
            datas.Add(new Data { ID = data.ID, X = x2, Y = y2, Y_general = Y2 });

            return datas;
        }

        private void ComputeL0AndN(double L, TransOption option, out double L0, out int n)
        {
            if (option == TransOption.Three)
            {
                n = (int)((L - 1.5) / 3) + 1;
                L0 = 3 * n;
            }
            else
            {
                n = (int)(L / 6) + 1;
                L0 = 6 * n - 3;
            }
        }

        private void ComputeAlphaCoeff(ref double[] alphaCoeff)
        {
            double e_square = Ellipsoid.e_square;
            double M0 = Ellipsoid.M0;

            alphaCoeff = new double[6];

            double A = 1 + 3 / 4 * e_square + 45 / 64 * Math.Pow(e_square, 2) + 175 / 256 * Math.Pow(e_square, 3) + 11025 / 16384 * Math.Pow(e_square, 4) + 43659 / 65536 * Math.Pow(e_square, 5);
            double B = 3 / 4 * e_square + 15 / 16 * Math.Pow(e_square, 2) + 525 / 512 * Math.Pow(e_square, 3) + 2205 / 2048 * Math.Pow(e_square, 4) + 72765 / 65536 * Math.Pow(e_square, 5);
            double C = 15 / 64 * Math.Pow(e_square, 2) + 105 / 256 * Math.Pow(e_square, 3) + 2205 / 4096 * Math.Pow(e_square, 4) + 10395 / 16384 * Math.Pow(e_square, 5);
            double D = 35 / 512 * Math.Pow(e_square, 3) + 315 / 2048 * Math.Pow(e_square, 4) + 31185 / 131072 * Math.Pow(e_square, 5);
            double E = 315 / 16384 * Math.Pow(e_square, 4) + 3465 / 65536 * Math.Pow(e_square, 5);
            double F = 693 / 131072 * Math.Pow(e_square, 5);

            alphaCoeff[0] = A * M0;
            alphaCoeff[1] = -1 / 2 * B * M0;
            alphaCoeff[2] = 1 / 4 * C * M0;
            alphaCoeff[3] = -1 / 6 * D * M0;
            alphaCoeff[4] = 1 / 8 * E * M0;
            alphaCoeff[5] = -1 / 10 * F * M0;
        }

        private void ComputeACoeff(double B, double X_zwx, ref double[] aCoeff)
        {
            double N = Ellipsoid.ComputeN(B);
            double t = Ellipsoid.ComputeT(B);
            double eta = Ellipsoid.ComputeEta(B);

            aCoeff = new double[7];
            aCoeff[0] = X_zwx;
            aCoeff[1] = N * Math.Cos(B);
            aCoeff[2] = 1 / 2 * N * Math.Pow(Math.Cos(B), 2) * t;
            aCoeff[3] = 1 / 6 * N * Math.Pow(Math.Cos(B), 3) * (1 - Math.Pow(t, 2) + Math.Pow(eta, 2));
            aCoeff[4] = 1 / 24 * N * Math.Pow(Math.Cos(B), 4) * (5 - Math.Pow(t, 2) + 9 * Math.Pow(eta, 2) + 4 * Math.Pow(eta, 4)) * t;
            aCoeff[5] = 1 / 120 * N * Math.Pow(Math.Cos(B), 5) * (5 - 18 * Math.Pow(t, 2) + Math.Pow(t, 4) + 14 * Math.Pow(eta, 2) - 58 * Math.Pow(eta, 2) * Math.Pow(t, 2));
            aCoeff[6] = 1 / 720 * N * Math.Pow(Math.Cos(B), 6) * (61 - 58 * Math.Pow(t, 2) + Math.Pow(t, 4) + 270 * Math.Pow(eta, 2) - 330 * Math.Pow(eta, 2) * Math.Pow(t, 2)) * t;
        }

        private void ComputeBCoeff(double Bf, ref double[] bCoeff)
        {
            double Mf = Ellipsoid.ComputeM(Bf);
            double Nf = Ellipsoid.ComputeN(Bf);
            double tf = Ellipsoid.ComputeT(Bf);
            double etaf = Ellipsoid.ComputeEta(Bf);

            bCoeff = new double[7];
            bCoeff[0] = Bf;
            bCoeff[1] = 1 / (Nf * Math.Cos(Bf));
            bCoeff[2] = -tf / (2 * Mf * Nf);
            bCoeff[3] = -(1 + 2 * Math.Pow(tf, 2) + Math.Pow(etaf, 2)) / (6 * Math.Pow(Nf, 2)) * bCoeff[1];
            bCoeff[4] = -(5 + 3 * Math.Pow(tf, 2) + Math.Pow(etaf, 2) - 9 * Math.Pow(etaf, 2) * Math.Pow(tf, 2)) / (12 * Math.Pow(Nf, 2)) * bCoeff[2];
            bCoeff[5] = -(5 + 28 * Math.Pow(tf, 2) + 24 * Math.Pow(tf, 4) + 6 * Math.Pow(etaf, 2) + 8 * Math.Pow(etaf, 2) * Math.Pow(tf, 2)) / (120 * Math.Pow(Nf, 4)) * bCoeff[1];
            bCoeff[6] = (61 + 90 * Math.Pow(tf, 2) + 45 * Math.Pow(tf, 4)) / (360 * Math.Pow(Nf, 4)) * bCoeff[2];
        }
    }
}
